      subroutine gmres(ncells, ijkcell, nvtx, ijkvtx, itdim, nsp, iwid, jwid, &
         kwid, precon, transpos, tiny, ibp1, jbp1, kbp1, c1x, c2x, c3x, c4x, &
         c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
         c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvolaxis, itsub, iter, error&
         , rnorm, fourpi, dt, cntr, clite, qdnc, qdnv, qom, bxv, byv, bzv, exmu&
         , eymu, ezmu, dive, residu, aq, q, bnorm, gradphix, gradphiy, gradphiz&
         , phi, diag, elle, nsmooth, wkl, wkr, wkf, wke, periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtx
      integer  :: itdim
      integer  :: nsp
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      integer , intent(in) :: itsub
      integer , intent(out) :: iter
      integer  :: nsmooth
      real(double)  :: transpos
      real(double)  :: tiny
      real(double) , intent(in) :: error
      real(double) , intent(out) :: rnorm
      real(double)  :: fourpi
      real(double)  :: dt
      real(double)  :: cntr
      real(double)  :: clite
      real(double) , intent(in) :: bnorm
      real(double)  :: elle
      real(double)  :: wkl
      real(double)  :: wkr
      real(double)  :: wkf
      real(double)  :: wke
      logical  :: precon
      logical  :: periodicx
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: vvol(*)
      real(double)  :: vvolaxis(*)
      real(double)  :: qdnc(*)
      real(double)  :: qdnv(itdim,*)
      real(double)  :: qom(*)
      real(double)  :: bxv(*)
      real(double)  :: byv(*)
      real(double)  :: bzv(*)
      real(double)  :: exmu(*)
      real(double)  :: eymu(*)
      real(double)  :: ezmu(*)
      real(double)  :: dive(*)
      real(double)  :: residu(*)
      real(double)  :: aq(*)
      real(double) , intent(inout) :: q(itdim,20)
      real(double)  :: gradphix(*)
      real(double)  :: gradphiy(*)
      real(double)  :: gradphiz(*)
      real(double)  :: phi(*)
      real(double) , intent(in) :: diag(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: nstart, n, ijk, ii, jj, k, iterp1
      real(double), dimension(21,2) :: g
      real(double), dimension(21) :: s
      real(double), dimension(21,21) :: h
      real(double) :: phibc, dtsq, dot, ws, g1, g2
!-----------------------------------------------
!
!
!
!
!
      phibc = 0.
      nstart = (ibp1 - 1)*(jbp1 - 1) + 1
      dtsq = dt**2
!
!     ZERO WORKING ARRAYS
!
      aq(ijkcell(:ncells)) = 0.0
!
      s(:itsub+1) = 0.0
      g(:itsub+1,1) = 0.0
      g(:itsub+1,2) = 0.0
      h(:itsub+1,:itsub+1) = 0.0
!
      q(ijkcell(:ncells),:itsub) = 0.0
!
!     CALCULATE THE INITIAL RESIDUAL ERROR
!
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phi, periodicx)
!
      call newres (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, jbp1&
         , tiny, cntr, nsp, itdim, dt, transpos, clite, qdnv, qom, bxv, byv, &
         bzv, exmu, eymu, ezmu, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, phi, gradphix, gradphiy, gradphiz, vvol, dive, qdnc, fourpi&
         , residu, elle)
!
      dot = dot_product(residu(ijkcell(:ncells)),residu(ijkcell(:ncells)))
!
      dot = sqrt(dot)
      s(1) = dot
!
      do n = 1, ncells
         ijk = ijkcell(n)
         q(ijk,1) = residu(ijk)/dot
      end do
!
!     ****************************************************************
!
!     begin gmres
!
!     ****************************************************************
!
      do iter = 1, itsub
!
!     apply preconditioner
!
         do n = 1, ncells
            aq(ijkcell(n)) = -q(ijkcell(n),iter)/diag(ijkcell(n))
         end do
!
!     compute A times preconditioned q
!
!dir$ ivdep
         do n = 1, ncells
            aq(ijkcell(n)) = phi(ijkcell(n)) + aq(ijkcell(n))
         end do
!
         call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, &
            phibc, aq, periodicx)
!
         call newres (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, &
            jbp1, tiny, cntr, nsp, itdim, dt, transpos, clite, qdnv, qom, bxv, &
            byv, bzv, exmu, eymu, ezmu, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x&
            , c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, &
            c6z, c7z, c8z, vol, aq, gradphix, gradphiy, gradphiz, vvol, dive, &
            qdnc, fourpi, aq, elle)
!
!
!dir$ ivdep
         do n = 1, ncells
            aq(ijkcell(n)) = residu(ijkcell(n)) - aq(ijkcell(n))
         end do
!
!
!    orthogonalize:
!
         do k = 1, iter
            dot = 0.0
            dot = sum(aq(ijkcell(:ncells))*q(ijkcell(:ncells),k))
            h(k,iter) = dot
!dir$ ivdep
            do n = 1, ncells
               aq(ijkcell(n)) = aq(ijkcell(n)) - dot*q(ijkcell(n),k)
            end do
         end do
!
!
         dot = sum(aq(ijkcell(:ncells))*aq(ijkcell(:ncells)))
         dot = sqrt(dot)
!
         iterp1 = iter + 1
         h(iterp1,iter) = dot
!
!     apply previous Givens rotations to h (update QR factorization)
!
         do k = 1, iter - 1
            ws = g(k,1)*h(k,iter) - g(k,2)*h(k+1,iter)
            h(k+1,iter) = g(k,2)*h(k,iter) + g(k,1)*h(k+1,iter)
            h(k,iter) = ws
         end do
!
!     compute next Givens rotation
         g1 = h(iter,iter)
         g2 = h(iterp1,iter)
         ws = sqrt(g1*g1 + g2*g2)
         g1 = g1/ws
         g2 = -g2/ws
         g(iter,1) = g1
         g(iter,2) = g2
!
!     apply g to h
         h(iter,iter) = g1*h(iter,iter) - g2*h(iterp1,iter)
         h(iterp1,iter) = 0.0
!
!     apply g to s
         ws = g1*s(iter) - g2*s(iterp1)
         s(iterp1) = g2*s(iter) + g1*s(iterp1)
         s(iter) = ws
!
!     |s(iter+1)| is the norm of the current residual
!     check for convergence
         rnorm = abs(s(iterp1))
         if (rnorm<=error*bnorm .or. iter==itsub) exit
!
!     normalize next q
!
         do n = 1, ncells
            q(ijkcell(n),iterp1) = aq(ijkcell(n))/dot
         end do
!
      end do
!
      s(iter) = s(iter)/h(iter,iter)
      do ii = iter - 1, 1, -1
         ws = 0.0
         ws = sum(h(ii,ii+1:iter)*s(ii+1:iter))
         s(ii) = (s(ii)-ws)/h(ii,ii)
      end do
!
!    compute new phi
      do n = 1, ncells
         ijk = ijkcell(n)
         phi(ijk) = -diag(ijk)*phi(ijk)
      end do
!
      do ii = 1, iter
!dir$ ivdep
         do n = 1, ncells
            ijk = ijkcell(n)
            phi(ijk) = phi(ijk) + s(ii)*q(ijk,ii)
         end do
      end do
!
      do n = 1, ncells
         ijk = ijkcell(n)
         phi(ijk) = -phi(ijk)/diag(ijk)
      end do
!
      return
      end subroutine gmres
